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ABSTRACT 

We describe our new " MLAPM-halo-finder" (MHF) which is based on the adaptive grid 
structure of the TV-body code MLAPM. We then extend the MHF code in order to track 
the orbital evolution of gravitationally bound objects through any given cosmological 
iV-body-simulation - our so-called "MLAPM- halo-tracker" (MHT). The mode of operation 
of MHT is demonstrated using a series of eight high-resolution A^-body simulations of 
galaxy clusters. Each of these halos hosts more than one million particles within their 
virial radii rvir- We use MHT as well as MHF to follow the temporal evolution of hun- 
dreds of individual satellites, and show that the radial distribution of these substruc- 
ture satellites follows a "universal" radial distribution irrespective of the host halo's 
environment and formation history. This in fact might pose another problem for simu- 
lations of CDM structure formation as there are recent findings by Taylor et al. (2003) 
that the Milky Way satellites are found preferentially closer to the galactic centre and 
simulations underestimate the amount of central substructure, respectively. Further, 
this universal substructure profile is anti-biased with respect to the underlying dark 
matter profile. Both the halo finder MHF and the halo tracker MHT will become part of 
the open source MLAPM distribution. 
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1 INTRODUCTION 

Over the last 30 years great progress has been made in the 
development of A'^-body codes that model the distribution of 
dissipationless dark matter. Algorithms have advanced con- 
siderably since the first A''^ particle-particle codes (Aarseth 
1963; Peebles 1970; Groth et al. 1977); we have seen the 
development of the tree-based gravity solvers (Barnes & 
Hut 1986), mesh-based solvers (Klypin & Shandarin 1983), 
then the two combined (Efstathiou et al. 1985) and multiple 
strands of adaptive and deforming grid codes (Villumsen 
1989; Suisalu & Saar 1995; Kravtsov, Klypin & Khokhlov 
1997; Bryan & Norman 1998; Knebe, Green & Binney 2001). 
While they all push the limits of efficiency in computational 
resources, each code has its individual advantages and limi- 
tations. The result of such research has been highly reliable, 
cost effective codes. However, producing the data is only one 
step in the process; the ensembles of millions of (dissipation- 
less) dark matter particles generated still require interpret- 
ing and then comparison to the real Universe. This necessi- 
tates access to analysis tools to map the phase-space which 
is being sampled by the particles onto "real" objects in the 
Universe; traditionally this has been accomplished through 
the use of "halo finders" . Halo finders mine Af-body data to 
find locally over-dense gravitationally bound systems, which 



are then attributed to the dark halos we currently believe 
surround galaxies. Such tools have lead to critical insights 
into our understanding of the origin and evolution of struc- 
ture and galaxies. To take advantage of sophisticated A'^- 
body codes and to optimise their predictive power one needs 
an equally sophisticated halo finder. 

Over the years, halo-finding algorithms have paralleled 
the development of their partner TV-body codes. We briefly 
outline the major halo finders currently in use: 

The Friends-of-Friends (FOF) (Davis et al. 1985; Frenk 
et al. 1988) algorithm uses spatial information to locate ha- 
los. Specifying a linking length bunk the finder links all pairs 
of particles with separation equal to or less than feiink and 
calls these pairs "friends". Halos are defined by groups of 
friends (friends-of-friends) that have at least one of these 
friendship connections. Two such advantages of this algo- 
rithm are its ease of interpretation and its avoidance of as- 
sumption concerning the halo shape. The greatest disadvan- 
tage is its simple choice of linking length which can lead to 
a connection of two separate objects via so-called linking 
"bridges". Moreover, as structure formation is hierarchical, 
each halo contains substructure and thus the need for dif- 
ferent linking lengths to identify "halos- within-halos" . There 
have been many variants to this scheme which attempt to 
overcome some of these limitations (Suto, Cen & Ostriker 
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1992; Suginohara & Suto 1992; van Kampen 1995; Okamoto 
& Habe 1999; Klypin et al. 1999). 

DENMAX (Bertschinger & Gelb 1991; Gelb & 
Bertschinger 1994a) and SKID (Weinberg, Hernquist & Katz 
1997) are similar methods in that thoy both calculate a den- 
sity field from the particle distribution, then gradually move 
the particles in the direction of the local density gradient 
ending with small groups of particles around each local den- 
sity maximum. The FOF method is then used to associate 
these small groups with individual halos. A further check is 
employed to ensure that the grouped particles are gravita- 
tionally bound. The two methods differ through their cal- 
culation of the density field. DENMAX uses a grid while 
SKID applies an adaptive smoothing kernel similar to that 
employed in Smoothed Particle Hydrodynamics techniques 
(Lucy 1977; Gingold & Monaghan 1977; Monaghan 1992). 
The effectiveness of these methods is limited by the method 
used to determine the density field (Gotz, Huchra & Bran- 
denberger 1998). 

A similar technique to the above is the Bound Density 
Maxima (BDM) method (Klypin & Holtzman 1997; Klypin 
et al. 1999) . In this scheme a smoothed density is derived by 
smearing out the particle distribution on a scale fsmooth of 
order the force resolution of the A''-body code used to gen- 
erate the data. Randomly placed "seed spheres" with radius 
^'smooth are then shifted to their local centre-of-mass in an 
iterative procedure until convergence is reached. Hence, as 
with DENMAX and SKID, this process finds local maxima 
in the density field. Bullock et al. (2001) further refined the 
BDM technique by first generating a set of possible centres, 
ranking the particles with respect to their local density and 
then implementing modifications which allow for credible 
identification of halos- within-halos. The Bullock et al. (2001) 
adaptation to BDM excels at finding halo substructure. 

When one is primarily concerned with distinct halos, 
all the mentioned methods perform exceedingly well. All ef- 
forts to refine and enhance those halo finding algorithms 
are due to the fact that A''-body codes overcame overmerg- 
ing only recently (Klypin et al. 1999) and are capable of 
finding satellites galaxies within dark matter host halos. It 
is therefore crucial to reliably identify "halos- within-halos" . 
In fact, one of the remaining problems for simulations of 
CDM structure formation is that high-resolution simulations 
nowadays predict far greater substructure (in total) than 
observed (Klypin et al. 1999; Moore et al. 1999). Results 
from gravitational microlcnsing suggest that the majority of 
substructure which does exist has to be close to the inner re- 
gions (Dalai & Kochanel 2002) which thus far has not been 
confirmed by such simulations. There are recent claims that 
although the ovcrmcrging problem has disappeared in the 
outer regions of the halo, the inner regions might still suffer 
from it (Taylor, Silk & Babul 2003). As these latter semi- 
analytic models do not suffer from such numerical problems, 
they find that such substructure docs exist in the inner re- 
gions. The question though arises as to whether there still 
remains an overmerging problem in the simulations or if cur- 
rent halo finding algorithms actually do break down at those 
scales. As we will discuss later, it becomes more difficult to 
locate peaks in the central region (if at all present) of the 
host halo due to a simple lack of contrast. 

In this paper we present a new method for identifying 
gravitationally bound objects in iV-body code output that 



uses the adaptive meshes of MLAPM (Knebe et al. 2001). This 
new code excelled at finding "halos-within-halos" revealing 
more substructure in the inner regions of the host halo. In 
its native form, our new algorithm works naturally "on-the- 
fly", but it has also been constructed with the flexibility 
necessary to handle a single temporal output from any N- 
body code. Our analysis software will become part of the 
publicly available MLAPM distribution*. The outline of the 
paper is as follows. In Section 2 we introduce the cosmo- 
logical models used to frame our discussion of the mode of 
operation of the new halo finder and tracker. A more detailed 
scientific analysis of this data set can be found in Paper II 
of this series (Gill et al. 2004a: hereafter, GKGDII). In Sec- 
tion 3 we introduce the new halo finder "MLAPM-halo-findcr" 
(MHF), describing its function, advantages, and limitations. 
Section 4 provides a brief analysis of the satellites found 
by MHF. In Section 5 we introduce the "MLAPM-halo-tracker" 
(MHT) which augments the halo finder by incorporating the 
ability to track the temporal evolution of satellites. Analysis 
of the halos tracked with MHT is described in Section 6. We 
next compaxe the two methods with other publicly available 
halo finding algorithms, such as FOF and SKID, in Sec- 
tion 7. We conclude with a summary and our conclusions in 
Section 8. 

This paper is the first in a series of three based upon 
the suite of simulations described herein. Paper II (GKGDII) 
investigates the satellite environments and their dynamical 
properties, while Paper III (Gill et al. 2004b) will investigate 
the tidal streams and debris from the disrupting satellites. 



2 SIMULATION DETAILS 

The A'^-body simulations presented in this and the compan- 
ion papers were carried out using the open source adap- 
tive mesh refinement code MLAPM (Knebe et al. 2001). MLAPM 
reaches high force resolution by refining high-density regions 
with an automated refinement algorithm. These adaptive 
meshes are recursive: refined regions can themselves be re- 
fined, each subsequent refinement having cells that are half 
the size of the cells in the previous level. This creates a hier- 
archy of refinement meshes of different resolutions covering 
regions of interest. The refinement is done cell-by-cell (in- 
dividual cells can be refined or de-refined) and meshes are 
not constrained to have a rectangular (or any other) shape. 
The criterion for (de-)refining a cell is simply the number of 
particles within that cell and a detailed study of the appro- 
priate choice for this number can be found elsewhere (Knebe 
et al. 2001). The code also uses multiple time steps on dif- 
ferent refinement levels where the time step for each level is 
a factor of two smaller than the time step on the previous 
level. The latest version of MLAPM also includes an adaptive 
time stepping that adjusts the actual time step after every 
major step to restrict particle movement across a cell to a 
particular fraction of the cell spacing, hence, improving the 
accuracy and computational time. 

We first created a set of four independent initial con- 
ditions at redshift = 45 in a standard ACDM cosmology 
(Qo = 0.3, = 0.7, V.bh'^ = 0.04, h = 0.7, as = 0.9). Next, 

* http : //astronomy . swin . edu . au/MLAPM/ 
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Table 1. Summary of the eight host dark matter halos. Distances are measured in /i ^ Mpc, velocities in km s ^, masses in 10^*/i ^ Mq, 
and the age in Gyrs. 



Halo 


Rvir 




Afvir 


^f'orm 


age 


A^sat(< r^ii) 


# 1 


1.34 


1125 


2.87 


1.16 


8.30 


158 


# 2 


1.06 


894 


1.42 


0.96 


7.55 


63 


#3 


1.08 


875 


1.48 


0.87 


7.16 


87 


#4 


0.98 


805 


1.10 


0.85 


7.07 


57 


#5 


1.35 


1119 


2.91 


0.65 


6.01 


175 


#6 


1.05 


833 


1.37 


0.65 


6.01 


85 


#7 


1.01 


800 


1.21 


0.43 


4.52 


59 


#8 


1.38 


1041 


3.08 


0.30 


3.42 


251 



512^ particles were placed in a box of side length 64/i~^ Mpc 

giving a mass resolution of rUp — 1.6 x 10^h~^ Mq. For each 
of these initial conditions we iteratively collapsed the closest 
eight particles to one particle reducing our particle number 
to 128^ particles. These lower mass resolution initial condi- 
tions were then evolved until z = 0. 

At « = 0, eight clusters from our simulation suite were 
selected in the mass range l-SxlO^^/i"^ Mq, each sampling 
differing environmental conditions. Then, as described by 
Tormen ct al. (1997), for each cluster the particles within 
two times the virial radius were tracked back to their La- 
grangian positions at the initial redshift (z = 45). Those 
particles were then regenerated to their original mass res- 
olution and positions, with the next layer of surrounding 
large particles regenerated only to one level (i.e. 8 times the 
original mass resolution), and the remaining particles were 
left 64 times more massive than the particles resident with 
the host cluster. This conservative criterion was selected in 
order to minimise contamination of the final high-resolution 
halos with massive particles. 

At the end of the high-resolution re-simulations the 

force resolution is dotorminod by the highest refinement level 
reached. The whole computational volume was covered by 
a regular domain grid consisting of 256^ cells. We had two 
separate criteria for refinement, a domain cell was refined 
when there was more than one particle per cell, further, ev- 
ery subsequent refinement was refined when there was more 
than four particles per cell. Thus the finest grid at z = 
consisted of 65,536 cells per side, giving a force resolution of 
w2/i,~^ kpc which allows us to resolve the host halos down 
to the central ~0.25% of the virial radii of the host halos 
(see Table 1). 

The halos chosen were selected to investigate the evo- 
lution of satellite galaxies and their debris in an unbiased 
sample of host halos, exploring the influence of environment 
upon the evolution of such systems. To achieve this goal, ex- 
cellent temporal resolution is required - as such we retained 
17 outputs from z — 2.5 to z = 0.5, equally spaced with 
At ^ 0.35Gyrs, supplemented with an additional 30 out- 
puts spanning z = 0.5 to « = with At « 0.17Gyrs. As we 
show in a companion paper (Gill et al. 2004), the average 
number of orbits for our satellites is of the order 1-2. There- 
fore we have approximately 10-20 outputs available to define 
the orbit of a satellite, which is more that adequate to follow 
a live orbit properly. We found that to sufficiently sample a 
live satellite orbit you need at least eight time-steps. As you 



increase the time sampling the stability of the result quickly 

converges. 

A simple analysis of the simulation at redshift z = 
provides us with the relevant information on the host halo. 
At 2 = the halo masses range from 1-3 xlO'^* Mq 
where the mass weis defined to be the total mass within 
the virial radius -Rvir, double counting both substructure 
and sub-substructure. The virial radii in turn were defined 
at the point where the mean averaged density of the host 
(measured in terms of the cosmological background density 
pb) drops below Avir = 340 with Mvir being the mass en- 
closed by that sphere. Wo then follow Lacoy & Cole (1994) 
and use their definition for formation time: the formation 
redshift Zform is the redshift where the halo contains half 
of its present day mass. Applying this criterion to our data 
we find that the ages of our host halos have a spread rang- 
ing from roughly 8.3 Gyrs to as young as 3.4 Gyrs. This 
alone shows that we are dealing with dynamically differ- 
ent systems even though their masses are comparable; our 
older halo's substructure has nearly twice the time to relax 
than the youngest one's satellites. A summary of the eight 
host halos is presented in Table 1 where the halos are pre- 
sented and numbered from oldest to youngest. The variation 
in the number of satellites from halo to halo with a trend 
for smaller hosts to contain less satellites can be accounted 
for by the mass cut applied to the satellites; as we expect 
a shift in the substructure mass function to lower masses 
for smaller hosts we are artificially cutting off satellites by 
applying a constant lower mass limit of 10^0h~^ Mg. 



3 MHF: MLAPM'S HALO FINDER 

The general goal of a halo finder is to identify gravitation- 
ally bound objects. As all halos are centered about local 
over-density peaks they are usually found simply by using 
the spatial information provided by the particle distribu- 
tion. Thus, the halos are located as peats in the density 
field of the simulation. To locate objects in this fashion, the 
halo finder is required in some way to reproduce the work 
of the Af-body code in the calculation of the density field or 
the location of its peaks. When locating halos like this, the 
major limitation will always be the appropriate reconstruc- 
tion of the density field. With that in mind we introduce 
MLAPM's-Halo-Findcr, MHF (or simply Finder) hereafter. 

MHF essentially uses the adaptive grids of MLAPM to lo- 
cate the satellites of the host halo. As previously mentioned 
in Section 2, MLAPM's adaptive refinement meshes follow the 
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density distribution by construction. Grid structure natu- 
rally "surrounds" the satellites, as the satellites are simply 
manifestations of over-densities within (and exterior) to the 
underlying host halo, a view which can best be appreciated 
through inspection of Figure 1. In this figure, the refine- 
ment grids of MLAPM are superimposed over the projected 
density of the particle distribution. The top image is the 5*'^ 
refinement level, with the 6"^ and 7*** levels shown below. 
We emphasise that the grids get successively smaller and 
are subsets of other grids on lower refinement levels. The 
advantage of reconstructing and using these grids to locate 
halos is that they naturally follow the density field with the 
exact accuracy of the A^-body code. No scaling length is 
required, in contrast with techniques such as FOF. There- 
fore, MHF avoids one of the major complications inherent to 
most halo finding schemes as a natural consequence of its 
construction. 

To locate appropriate halos within our simulation out- 
puts we first build a list of "potential centres" for the halos. 
Using the full adaptive grid structure invoked by MLAPM, with 
the same refinement criterion as for the original runs, we re- 
structure the hierarchy of nested isolated MLAPM grids into a 
"grid tree" and generate a list of prospective halo centres by 
storing the centroid of the densest grid at the end of each 
grid tree's "branch". Assuming that each of these peaks in 
MLAPM's adaptive grids is the centre of a halo, we step out in 
(logarithmically spaced) radial bins until the density reaches 
Psatoiiitc(r-vir) = A vir (2)pb (z) , where pb is the universal back- 
ground density, unless we reach a point rtrunc where an up- 
turn in the radial density profile is detected. This rise is 
encountered for (almost) all satellites embedded within the 
background density of the host halo, a point that we will 
discuss in more detail in Section 4. The outer radius of the 
satellite is defined to be either rvir or rtrunc, whichever is 
smaller, and dubbed rme- Using all particles interior to rme 
we calculate other canonical properties for each halo such as 
its mass, rotation curve, and velocity dispersion. 

We now, however, need to prune the list of (still prospec- 
tive) halos by removing gravitationally unbound particles 
and duplicate halos. The latter occurs in two steps - first, for 
each satellite a set of "duplicate candidates" is constructed 
based on the criterion that their centres lie within each oth- 
ers' outer radii r-me- Second, this list is then checked by com- 
paring the internal properties of the candidates. A candidate 
was affirmed to be a duplicate once its mass, velocity dis- 
persion, and center of mass velocity vector agreed to within 
80%. We then kept the halo with the higher central density 
and removed the other one from the satellite catalogue com- 
pletely. This is a rare circumstance, yet one to which we will 
return in Section 4. With our nearly complete set of halos 
now in hand, we proceed to remove gravitationally unbound 
particles. This again is done in an iterative process. Starting 
with the MHF halo centre, we calculate the kinetic and po- 
tential energy for each individual particle in the respective 
reference frame and all particles faster than two times the 
escape velocity are removed from the halo. We then recalcu- 
late the centre, and proceed through the process again. This 
pruning is halted when a given halo holds fewer than eight 
particles or when no further particles need to be removed. 
We finish by recalculating the internal properties of the ha- 
los with the radial density profiles of the satellites fitted to 



t « 






Figure 1. This panel shows a series of 3 consecutive refinement 
levels of MLAPM's grid structure starting at the 5th refinement 
level superimposed upon the density projection of the particle 
distribution. 



the functional form proposed by Navarro, Frenk & White 
(1997; hereafter, NFW) 



M{<r) 



(1) 



(r/r,)(l + r-/r,)2' 
in the range from 8h~^ kpc (~ 4 x force resolution) to VmF 
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The scale radius Vs is used to define the concentration of the 
halo 

c = rvir/rs. (2) 

The procedure outlined above naturally deals with over- 
lapping halos and substructure halos, respectively. But as 
mentioned before, for such objects the virial radius can not 
be determined properly as we will observe a rise in the radi- 
ally binned density profile due to the overlap with another 
halo or the embedding into the host. In that case we set the 
outer radius of the (sub-)halo to be that point where the 
density profile rises and all canonical properties are derived 
using all (gravitationally bound!) particles interior to that 
radius. And the fit to an NFW profile Equation (1) is only 
done out to that radius, too. The situation is different once 
both of the overlapping halo's centres are within each other's 
virial/upturn radius: we then checked, if those two objects 
are just duplicates by comparing their internal properties. 

As stated in Section 2, it is our aim to investigate the 
evolution of satellite galaxies within their host halos. Thus, 
we restrict our satellites to having at least 50 high-resolution 
simulation particles, which corresponds to a mass-cut of 
Afcut ~ W^'^h~^ M©. Moreover, each satellite must contain 
at least 50% percent of its mass in high-resolution particles. 
In practice, this latter constraint is not a critical one, rel- 
evant only for satellites beyond twice the host halo's virial 
radius. 

We can further take advantage of the MLAPM grids for 
measuring the triaxiality of regions surrounded by an iso- 
density contour. Essentially the various refinement levels are 
cuts in the density field (isodensity surfaces) . We calculated 
the inertia tensor for each isolated refinement, weighting 
each cell by its density. Then using the eigenvalues of the 
inertia tensor we construct the triaxiality parameter (Franx, 
lUingworth & Zeeuw 1991) 

T = (a" - b^)/{a' - c") . (3) 

To describe the host halo's triaxiality we used the 6"^ 
refinement level in MLAPM. According to the refinement crite- 
rion adopted in the simulations the 6**^ level surrounds ma- 
terial about 3000 times denser than pb or, in other words, 
nine times denser than the material at the virial radius. A 
density of roughly 9 x p(rvir) corresponds to approximately 
the half-mass radius of the host. 

MHF is implemented into MLAPM in a way that provides 
the user simultaneously with a snapshot of the dark matter 
particles and halo catalogues at each required output. The 
most obvious advantages of having the analysis performed 
"on the fly" are the reduction in computer and human hours 
in the initial halo analysis stage. Embedding the halo anal- 
ysis in the code also enables us to potentially analyse the 
data at unprecedented time resolution, if required. 1^ How- 
ever, MHF can also be used with any already existing single 
time-step snapshot and hence is not limited to data pro- 
duced by MLAPM; it can also be used for any A''-body output 
provided the latter is converted to MLAPM's binary format 
using the tools included in the MLAPM distribution. 



t MHF can be switched on either to act only when writing an 
output file (-DMHF) or at each individual time step (-DMHFstep). 
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Figure 2. Number of satellites (normalised) orbiting within the 
virial radius of the host halo at 2 = 0, as a function of radial 
distance r, normalised by the virial radius of the respective host 

4 ANALYSIS OF MHF HALOS 

MHF was applied to each of the 376 temporal outputs (47 
outputs per each of the eight independent halos), providing 
us with a list of all satellites and their internal properties 
at each individual redshift under consideration. As stated 
earlier, the detailed analysis of the science associated with 
this study is presented in Paper II (GKGDII). We do however 
wish to highlight several key preliminary results here which 
relate specifically to the halo identification process. 

In Figure 2 we plot the normalised number of satel- 
lites as a function of normalised radius. One, perhaps not 
surprising, aspect of Figure 2 is the similarity in the slopes. 
Although the number of satellites in each halo may vary, the 
relative radial distribution of the satellites is similar across 
halos. This is reminiscent of the universal density profile 
of dark matter halos, as described by NFW. Although the 
radial distribution of the satellites remains consistent, there 
exists a range of substructure densities for the halos, as there 
is a spread in the number of satellites within each halo (re- 
call Table 1). Therefore, we should be able to distinguish the 
effects of substructure density on the physical properties of 
the satellites. 

The other striking feature of Figure 2 is the lack of satel- 
lites in the inner 15% of the virial radius. One might ask why 
this is the case. Does it indicate that satellites dynamically 
avoid the central region of the halo? Perhaps they simply 
do not spend much time there? Perhaps via physical means 
satellites that venture near the centre are either merged or 
experience such strong tidal forces that they are destroyed? 
Perhaps we are simply dominated by numerical effects and 
are witnessing the premature destruction of halos in dense 
environments. This latter problem - known as overmerging 
- affected low-resolution dissipationless simulations in the 
early 1990s, failing to produce galaxy-sized dark matter ha- 
los in clusters (e.g. Summers et al. 1995; van Kampen 1995; 
Moore et al. 1996). Traditionally, this was explained by the 
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lack of dissipation in the simulations. With the inclusion of 
a baryonic component, denser objects could form and sur- 
vive in the centres of these dense regions. However, with the 
onset of higher resolution simulations a converse effect was 
encountered - specifically, an abundance of substructure was 
found (Klypin et al. 1999). 

The explanation of overmerging (or substructure dis- 
ruption) was accredited to numerical limitations in the sim- 
ulations, van Kampen (1995) found that particle evapora- 
tion due to two-body effects is only important for low par- 
ticle number halos (<30 particles ). Moore et al. (1996) 
further investigated particle halo heating, which they con- 
cluded was negligible should sufficient mass resolution ex- 
ist. Moore et al. also demonstrated the for satellites in a 
static host potential, if a simulation had insufficient spa- 
tial resolution, halos would have artificially large cores and 
hence undergo accelerated tidal disruption. They also found 
that halos become unstable and arc erased when the tidal 
radius is smaller than approx 2-3 times the halo core ra- 
dius (which itself can be related to the gravitational soft- 
ening length). Klypin et al. (1999) investigated the issue of 
"overmerging" in great detail using a variety of higher res- 
olution simulations, concluding that the resolution required 
to avoid artificial destruction of galaxy-sized halos of mass 
10"/i"^ Mo was < 2 kpc (spatial) and < lO^/i"^ M© 
(mass) . 

Since we appear to have sufficient numerical resolution 
and our data lies well within the limits of not being domi- 
nated by overmerging, one might query whether or not the 
lack of substructure in the inner region is due to a limi- 
tation of our halo finder. When defining the radius of our 
halos we could not for all halos follow the density profiles 
out to Tvir defined via Psat(?'vir) = AvirPi,; as noted earlier, 
it was necessary in many cases to define a truncation ra- 
dius rtrunc. The existence of nrunc generally indicates that 
the satellite is embedded within the host's density field, as 
already noted by Bullock et al. (2001). Thus, as a satellite 
gets closer to the central density region of the host halo, its 
overdensity peak becomes less contrasted. It is intrinsically 
harder to find satellites with low central densities under the 
standard paradigm of halo finding, especially close to the 
cuspy centre. It is not at all obvious how to disentangle the 
particle distribution of the satellite and the host halo: this is 
a fundamental limit to finding halos in the traditional way of 
observing over-densities and requires further investigation. 
In the next section, however, wo introduce a method of find- 
ing halos that eliminates the background halo and, hence, 
minimises this problem. 

The MHF method fails in the inner regions for two rea- 
sons. Firstly, because it is hard to detect the upturn in the 
density field, substructure is eliminated through suspected 
duplication of a halo because the substructure's radius has 
been falsely tracked out to essentially the virial radius of the 
host, its own upturn radius has been missed. The second 
reason results from a fundamental flaw in MHF's methodol- 
ogy - that the smaller satellite grids merge with the host's 
refinement grid and hence do not produce an isolated refine- 
ment. Therefore we are losing potential centres, a problem 
illustrated further in Figure 3. There, we show the inner 
250/i"^ kpc of halo #1, along with the grids for the 7"* re- 
finement level of MLAPM (gray-shaded areas) with the central 
refinement about 60/i~^ kpc in radius. The dark spheres in- 



dicate the positions of satellites located by MHF. Note that 
those dark spheres that do not encompass an isolated re- 
finement grid would do so at one of the next coarser (or 
finer) levels. The light sphere at the border of the ellipsoidal 
host refinement also surrounds a satellite galaxy. However, 
this object was not picked up by the Finder but rather by 
the Tracker outlined in the next section. MHF was unable to 
identify this satellite as an individual object as its refinement 
grid has merged with the host's grid, thus not allowing an 
isolated refinement and a potential center, respectively. The 
straight line pointing to this satellite is simply its orbital 
path. 

The problem can be viewed differently in Figure 4. Here 
we plot the radially averaged density of the host halo at the 
position of a satellite against the maximum, central density 
of the satellite itself for redshift z = Q. The results are pre- 
sented for the Finder (crosses) as well as for the Tracker 
(diamonds) to be introduced in Section 5. The line running 
through the plot for each individual host halo marks the 
1:1 correspondence: satellites that fall onto (or even above) 
this line have central densities equal to (or smaller than) the 
host environment they are embedded within. We do observe 
a general (and reasonable) trend for satellites to have higher 
central densities than their dark matter vicinity. However, 
the figure also proves that the Tracker tends to also find 
satellites less contrasted and closer to the 1:1 relation, re- 
spectively. When interpreting Figure 4, and especially com- 
paring the Finder to the Tracker results, one needs to bear 
two things in mind: firstly, there are many more Tracker- 
satellites obscuring a one-to-one comparison with Finder, 
and secondly, Finder relies on rme as the final point of the 
profile whereas Tracker has the ability to properly measure 
Tvir. Minor differences in binning can also lead to very small 
changes in the central density calculation. 



5 MHT: MLAPM'S HALO TRACKER 

Conventional halo finders have a rich history in identify- 
ing isolated systems. In this regard, MHF might be viewed 
as simply an alternative approach to an already reasonably 
well-understood problem. To be fair though, MHF does push 
the conventional paradigm of simply using the three dimen- 
sional spatial data to locate the halos to the limit, locating 
the halos with (nearly) the exact accuracy of the iV-body 
code. Having this ability, MHF becomes the ideal halo finder 
to locate substructure for MLAPM and no doubt an excellent 
halo finder for other codes. Although, as we have seen in 
Section 4, apart from numerical limitations in A^-body codes 
(e.g. overmerging) there still remain limitations to our cur- 
rent halo finders. These limitations become a problem in the 
simulations when considering the substructure of any dense 
system, for example galaxy clusters and galaxies. Therefore, 
to successfully find the substructure we need to change the 
paradigm used to find it. 

To successfully make this change we must first under- 
stand the environment in which we are finding the substruc- 
ture and then exploit its characteristics. One characteristic 
is that most halos conserve their identities, that is substruc- 
ture halos rarely undergo mergers in halos because of their 
high relative velocities (Ghigna et al. 1998; Okamoto & Habe 
1999). Further, halos in dense environments undergo tidal 



© 0000 RAS, MNRAS 000, 000-000 



Identifying dark matter substructure 7 




10' 



10 



Figure 3. The inner 250 Icpc of halo #1 with the particles' 
line of sight density shown. Wc show the grids of the 7*'' refine- 
ment level of MLAPM, with the central refinement about 60h^^ kpc 
in radius. The dark spheres represent the satellites located by MHF. 
The light sphere surrounds a satellite galaxy not found by MHF. 
The apparent sizes of the spheres are simply a visualisation eft'ect 
as spheres farther away from the virtual observer appear smaller. 



stripping and substructure interactions, no longer accreting 
material, but being stripped of it. Thus in such environments 
it is sufficient to trace the particles of the satellite once the 
satellite has entered the host's virial radius. In this section 
wc introduce MLAPM's-Halo- Tracer, MHT (or simply Tracker) 
hereafter. 

MHT takes an arbitrary output of our new MLAPM-based 
halo finder MHF, and correlates the particles for an arbitrary 
number of time steps using all simulation outputs from that 
initial output until rodshift z — 0. In our particular inves- 
tigation it was appropriate to define that initial arbitrary 
output to be the formation time Zfoim of the host halo or 
the time when the host halo contained half of its present 
mass (Lacey & Cole 1994). We then followed all the satellites 
that were within two times the virial radius of the host halo 
at this formation time. Although, we miss a few satellites 
due to MHF's identification limitations in the inner 10-15% 
of the halo, from this time on MHT precisely follows the or- 
bits of our initial sot of satellites irrespective how close they 
come to the host's centre. Explicitly, MHT takes the particles 
from an initial MHF analysis and then locates these parti- 
cles in the next available output again. Tracker's first task 
is then to (rc-)calculatc the halo's centre. This is done by 
using the centre-of-mass of the innermost 20 particles from 
the previous time step as an initial estimate, then using the 
same iterative method to check the credibility of the halo, 
as outlined in Section 3. Once the satellite was identified 
as bona fide the radial profile was generated, a NFW pro- 
file fitted, and other canonical properties calculated. The 
binning for the profile again used logarithmically spaced ra^ 
dial bins covering the entire particle distribution. The ra^- 
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Figure 4. The density of the host halos at the radial distance 
-Dsat of the satellite versus the central density of the satellites. 
The crosses represent the satellites found by Finder while the 
diamonds are the satellites found by Tracker. 



dius of the satellite was consistently determined as being 
the radius when the cumulative density profile dropped be- 
low psatciiitc(?'vir) = ^\-\i{z) Pb{z) ■ Tliis timo wc wiU not en- 
counter the situation where the profile rises again as is the 
case for Finder; the satellite is no longer embedded within 
the "particle background" of the host halo but treated as a 
separated entity. 

There arc a number of advantages in tracing the halos 
in this way - first, because we are tracking just the satellites' 
individual particles we do not have the complication of the 
background density distribution and the consequent lack of 
contrast against the host system for all outputs z < 2:form. 
Following from this, we do not have to accept the truncation 
radius - the radius where the Finder encounters and upturn 
in the denisty profile - as the "natural" radius of the satellit. 
Further, this method allows us to investigate the develop- 
ment of tidal streams, which forms the basis of the extensive 
analysis provided in Paper III. 



6 ANALYSIS OF MHT HALOS 

The nature of hierarchical structure formation, i.e. mergers, 
dynamical and tidal destruction of substructure, requires a 
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little more work when applying MHT to our simulation data. 
Even though were are now tracing the initially bound parti- 
cles forward in time, we need a criterion to decide whether 
a satellite galaxy is disrupted or still alive. We therefore in- 
troduce the tidal radius as given below. 

6.1 Tidal Radius 

The tidal radius is defined to be the radius of the satellite 
where the gravitational effects of the host haJo are greater 
than the self-gravity of the satellite. When approximating 
the host halo and the satellite as point masses and main- 
taining that the mean density within the satellite has to be 
three times the moan density of the host halo within the 
satellites distance D to the host halo (Jacobi limit) the def- 
inition for tidal radius reads as follows 



nidai 



D 



(4) 



V3m; 

where m is the mass of the satellite and M = M(< D) is 
the mass of the host halo internal to the distance D. 

In order to stabilise the determination of the tidal ra- 
dius rtidai we actually use an iterative procedure again. By 
defining the satellite mass m as the mass internal to rtidai, 
i.e. m = m{< nidai), we find nidai by solving 

3 

D = . (5) 



nidal 



m{< rtidai) 



3M 



Starting with the distance of the furthest particle in the 
satellite particle distribution as the initial guess for rtidai, the 
method quickly converges in only two-to-three iterations. 

6.2 Satellite Disruption 

As the satellites orbit within the host halo they undergo 
tidal stripping, hence, the satellites arc gradually losing 
mciss. Therefore, their particle distribution becomes more 
and more diffuse, reducing the tidal radius of the satellite. 
Eventually there comes a point when the satellite loses suffi- 
cient mass that we begin to reach the limits of our numerical 
resolution, not having sufficient rmmber of particles to fol- 
low the satellite further. The satellite might have survived 
for longer, however, we do not have the resolution to fol- 
low it. Thus when a satellite passes through our numerical 
limit we tag it "disrupted". In practice this means that if 
there are fewer than 15 particles within the tidal radius we 
classify the satellite as being disrupted. Note that we are 
unable to separate numerical resolution disruption and real 
physical disruption of a satellite. It is not clear if we had in- 
finite mass resolution that the satellite would still actually 
survive. 

We also stress that since our tidal radius formulae as- 
sumes circular orbits, the usual D corresponds to the peri- 
centric distance (Hayashi ot al. 2003); however, wo calculate 
rtidai for each satellite at each individual output to check for 
(tidal) disruption. 

Furthermore, particles outside the tidal radius are not 
automatically stripped from the halo - what is just as im- 
portant is the time spent under the influence of that tidal 
field, it still takes time for the particle to climb out of the 
potential. For example a satellite might be on a very eccen- 
tric orbit and pass close to the centre of the host halo, thus 





0.4 0.6 
redshift z 



Figure 5. Distance from the center of the host halo in real co- 
ordinate system as a function of redshift for one particular satel- 
lite with (Msat/Mhost = 0.7 X IQ-^). 



having a very small tidal radius at the pericentre. Now it is 
true that most of the tidal stripping will occur at this peri- 
centre, however, because the satellite only spends a short 
time there not all the particles outside the tidal radius will 
be stripped. 



6.3 Orbital Information 

In GKGDII we present a detailed analysis of the satellites' 
orbits, but do take the opportunity here to provide some ba- 
sic terminology relevant to both Papers I and II here. The 
high temporal resolution of our simulations {At — 0.17 and 
At = 0.35 Gyrs, respectively, for z > 0.5 and z < 0.5) en- 
ables us to track in detail the orbits of the satellites. As an 
example, in Figure 5 we show the orbit of one particular 
satellite. We can see that this satellite initially plunged in 
from outside the virial radius a,t z — 0.8 and was subse- 
quently captured by the host, undergoing two further orbits 
prior to z = 0. This orbital information was then used to 
construct a measure of eccentricity 



1-P 



(6) 



where p is the most recent "closest" distance to the host's 
centre as a minima (labelled peri in Figure 5) and a the 
most recent "furthest" distance (labelled apo) as a maxima. 
Moreover, we are also in the position to calculate the num- 
ber of orbits from the time evolution of the distance to the 
host centre. To this extent, we simply count the number of 
extrema (four in the case shown) and divide by two. We fur- 
ther correct for incomplete orbits at the beginning and end 
points of the distance relation, resulting in the following for- 
mula for the total number orbits: 



Nr.- 



Nr. 



:trcma . . / 1 ^1 \ , • / ^ \ /„\ 

-^min(-,-)-Fmin(-,- ) (7) 

Z Z 12 ^ tn — 1 
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Figure 6. Number of satellites within a particular radius nor- 
malised by the total number of satellites as a function of radius 
normalised by the virial radius at 2 = 0. Only satellites more 
massive than 2 X 10^^h~^ Mq were taken into account. 



Figure 7. Number density of satellites normalised by the average 
number of satellites per unit volume as a function of radius nor- 
malised by the virial radius at z = 0. Only satellites more massive 
than 2 X 10^^h~^ M0were taken into account. 



The number of orbits measured by that method for the sam- 
ple satellite presented in Figure 5 is A'orbits = 2.69. 

We emphasise though that the orbits of the satellites 
are not always as aesthetically "smooth" as that for the one 
presented in Figure 5. We are dealing here with live poten- 
tials and hundreds of satellites orbiting within it simulta- 
neously. The host halo is constantly growing in mass and 
shows internal oscillations in shape due to ongoing mergers 
(see GKGDII). This has has an impact on the orbital evolu- 
tion of the satellites, as described in Paper II. 



6.4 The radial distribution of satellites 

In Figure 2 we showed the radial distribution of satellites 
for MHF. We now present in Figure 6 the same plot for MHT 
highlighting the superiority of the Tracker. Once again we 
observe the similarity in the slopes across the eight halos, a 
so-called "universal satellite distribution". However, the im- 
portant result is that using the Tracker we now find (more) 
objects within the central 10% of r^ii of the host halos. 
For halo #1, for instance, we located 5 satellites with mass 
greater that 10^°/i~^ Mq within 10% of the virial radius 
with the closest satellite at 2 = being a mere 35/1"^ kpc 
away from the host's centre. 

To allow a more natural comparison to work published 
by other authors, we present the data in a slightly different 
fashion in Figure 7. This time the radial number density of 
satellite galaxies is shown. As with Ghigna et al. (2000) and 
De Lucia et al. (2003) we also find that the subhalo popula- 
tion is "anti-biased" relative to the dark matter distribution 
in the inner regions of the halos. Moreover, we again observe 
no trend with environment for the sample of eight halos un- 
der investigation; all halos, irrespective of age and richness, 
do show the same anti-bias in the satellite distribution. 



7 COMPARISON TO OTHER HALO FINDERS 

In this section we compare MHF and MHT to two other halo 
finders, namely SKID and FOF. In Figure 8 we present a 
visual comparison of the effectiveness of the respective halo 
finders. Firstly, the top two panels show the line of sight den- 
sity projection of particles for halo #1 within a 700h~^ kpc 
sphere. The left panel displays all the particles, while the 
right panel only shows every third particle. In the remaining 
four panels, we present the results of the halo finding algo- 
rithms: (reading clockwise, starting upper left) MHT, SKID, 
FOF, and MHF. A sphere of fixed radius 20h~^ kpc surrounds 
each located satellite where different apparent sizes are sim- 
ply visualisation effects, i.e. spheres farther away from the 
virtual observer appear smaller. 

SKID was run with multiple linking lengths and the op- 
tion to remove gravitationally unbound particles enforced. 
Under close visual inspection, the best SKID results were 
found with linking lengths 6=0.03 & 6=0.05 times the inter- 
particle separation. The same values were taken for the FOF 
analysis. The results from the analysis with linking length 
6=0.05-^ appear to be the most reliable ones and hence are 
shown in Figure 8. We need to stress though that the anal- 
ysis of the SKID and FOF results can be further refined by 
combining various halo catalogues into a tree as explained 
in, for instance, Ghigna et al. (2000). However, one of the 
benefits our MLAPM halo finders is that neither requires any 
further input such as a linking length. Moreover, the lat- 
est version of MLAPM is capable of performing the analysis 
"on-the-fly" . A visual inspection of Figure 8 indicates that 

t For a detailed description of SKID please refer to 
http : //www-hpcc . astro .Washington. edu/tools /skid. html. 
§ A linking length of 6=0.05 expressed in terms of the inter- 
particle separation translates into a physical linking length of 
r^6h~^ kpc, which is roughly three times the force resolution. 
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Figure 8. The top two panels show the hne of sight density 
projection of particles for halo #1. The left panel displays all 
the particles, while the right panel shows every third particle. 
The remaining four panels show the results of the halo finding 
algorithms: MHT, SKID, FOF, and MHF (clockwise from top left). 
A sphere of fixed radius 20/i~^ kpc surrounds each satellite. 

MHT provides the most complete halo hst. Specifically, the 
Tracker found 53 satellites, the Finder found 32, SKID 33 
and FOF 32 within the plotted spherical region of diameter 
lAh~^ Mpc. Within the sets of satellites, there is of course 
considerable overlap. Essentially, each set of satellites found 
by MHF SKID, and FOF were subsets of MHT. To look at 
this quantitatively, we calculate the radial number density 
of satellites for all four halo finders again (cf. Figure 7). This 
time we concentrate on halo #1 highlighting the differences 
between the halo finding methods. The result is presented in 
Figure 9 which clearly shows the success of the Tracker over 
all other methods. However, it is interesting to note that a 
simple FOF analysis gives quantitatively similar results to 
the more sophisticated SKID data. The difference between 
the Finder and the Tracker is quite remarkable, as is the 
similarity between MHF and SKID with 6=0.03. Our explana- 
tion for the lack of substructure in the central region for this 
particular SKID analysis is that these objects where either 
removed because of our lower mass cut of 100 particles or 
the fact that SKID did not classify them as gravitationally 
bound. What we can also learn from Figure 9 is that the sen- 
sitivity of the halo finder in the inner regions can severely 




0.1 1.0 



Figure 9. Number density of satellites normalised by the average 
number of satellites per unit volume as a function of radius nor- 
malised by the virial radius at z = 0, halo #1 for each of the halo 
finding algorithms: MHT, SKID, MHF, FOF. Only satellites more 
massive than 2 X lQ^^^h~^ MQwere taken into account. 
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Figure 10. Same as Figure 4, but this time comparing Tracker 
(diamonds) and Finder (crosses) with the SKID, b=0.05 (trian- 
gles) analysis for halo #1. 

bias the results. For example using MHF or SKID with b=0.03 
provides a much stronger anti-bias in the satellite distribu- 
tion than the more appropriate Tracker and SKID (6=0.05) 
analysis, respectively. 

We like to the close the comparison by coming re- 
examining Figure 4, this time including the results from 
the SKID (6=0.05) analysis. The result for halo #1 can be 
viewed in Figure 10. We still observe that only the Tracker 
is capable of resolving satellites with central densities close 
to the (local) density of the host halo. We also inspected the 
situation for the FOF (6=0.05) data and could not find any 
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significant difference and hence decided to not plot the data 
for clarity. 



8 CONCLUSIONS 

Computational cosmology is not only limited by crucial fac- 
tors such as the dynamical range and the mass resolution, 
but also by its analysis tools. We emphasise, perhaps obvi- 
ously to most readers, that AT-body codes simulating struc- 
ture formation in the Universe can only ever be as useful as 
their associated analysis tools allow. 

In this paper wc presented two new methods for iden- 
tifying gravitationally bound objects in such simulations. 
Both methods are based upon the open source adaptive 
mesh refinement code MLAPM (Knebe et al. 2001). They 
both exploit the refinement hierarchy of said code and 
hence locate halos as well as halos-within-halos with exactly 
the same accuracy as MLAPM simulates their evolution. We 
showed the limitations of a simple snapshot analysis and 
how it can be overcome by taking into account the whole 
history of each halo. Thus not restricting the halo finding 
algorithm to just the spatial information but rather includ- 
ing the velocity information as well, hence, using the full 
six-dimensional information available. 

Not only do we intend to implement the halo finder into 
the distribution of MLAPM allowing for an on-the-fly analy- 
sis saving both computational and human resources in the 
analysis process, but it is also our intention to have it as 
a stand alone program. In both cases, the implementation 
will be such that it can analyse a single output of any given 
A'^-body code. 

We showed that halo finding still possesses the inherent 
problem of overmerging in the very central regions of the 
host. However, by tracking satellites rather than identify- 
ing them at separate time snapshots of the simulation wc 
learned that this problem is not overmerging in the conven- 
tional sense. The objects are in fact present and simulated 
properly, but their densities have insufficient contrast to be 
picked up by a simple "finder" . Only when tracking them in 
time from (at best their very own) formation time were we 
able to quantify their existence as close as 5% of r^°/* to the 
cuspy centre. We do not intend to question the credibility 
of other (most excellent) halo finders such as SKID though. 
We rather pointed out that the results in the central region 
are subject to subtleties that can be most easily avoided by 
tracking satellites. 

It has recently been pointed out by Taylor et al. (2003) 
that the radial distribution of the Milky Way's satel- 
lite galaxies does not reconcile with the predictions of 
semi-analytical models of galaxy formation and cosmolog- 
ical simulations, respectively, which has been confirmed by 
Kravtsov et al. (2004) . Even though there is quite a promi- 
nent substructure population in the simulations it is mostly 
clustered in the outer regions of the host whereas the Milky 
Way satellites are preferentially found closer in. It remains 
unclear if this poses a new challenge to the CDM structure 
formation scenario or simply reflects what we presented in 
this study. Namely, identifying substructure halos that lie 
close to the centre of the host is intrinsically challenging 
and might be overcome by actually tracking the satellites 
rather than finding them. 



Finally, we further increase the statistics to suggest that 
anti-bias of the satellite distribution is a common property of 
the satellite distribution, perhaps even universal. However, 
this antibias is still most likely a relic of numerical over- 
merging. A convergence study for the satellite populations 
is critical in resolving this issue and reconciling the lens- 
ing observations. However, one very interesting result is the 
continued self-similarity in the dark matter distribution, the 
satellite radial distribution is common, once again perhaps 
"universal" even though the number of substructure satel- 
lites changes, reminiscent of the universal density profile of 
the underlying host. 
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